doeverything= function() {


 
totalVar = 0.4899106
totalSD = sqrt(totalVar)
propGroupVar = 0.8


parameters = list(
	   geno=log(1), env=log(1), "geno:env"=log(1),
		 size=150000, probGeno=0.2, 
		 probEnv=0.3, genoMiss= 0.1,
		 envMiss = 0.1, oneEnvPerCommunity =T,
     cancerMissRate=0.15, falseCancerRate=0.00045, 
		 sdGroup= sqrt(propGroupVar) * totalSD, 
		 sdPerson=0#,sqrt(1-propGroupVar)*totalSD, 
		 Enrollmentprobs = c(20000, 40000, 50000, 40000), 
		 agerange = c(35, 69),
		 Followup = c(10, 20, 30),
		 Ncommunity =  50,
		 EqualCommunity=T
)


#populationData = getPopData("../data/On129Communities.csv")  
populationData = getPopData()  


deathData = list(
  "F"=read.table("../data/MortalityNocancerFemale.txt", header=T), 
  "M"=read.table("../data/MortalityNocancerMale.txt", header=T)
)
DiffCancer= list(
	"F"=lambdaDiffCancer(myfile="../data/IncidentFemaleDiffCancer.csv",
	 	agegroup),
   M=lambdaDiffCancer(myfile="../data/IncidentMaleDiffCancer.csv",
	 	agegroup)
	 )	
  

lambda =list("F"=list(x=seq(30, 90, by=5), y= exp((-0.5)*0.585)*c(116.4,169.7,272.6,401.8,555.4,746.8,1003.2,1254.8,1528.7,1741.4,
1903.1,1950.2,1922.5) / 100000)  ,
M=list(x=seq(30,90, by=5), y=exp((-0.5)*0.585)*c(61.8,81.2,123.6,247.8,466.8,877.5,1403.7,2096.6,2598.4,2887.2,
3086.4, 3156.3, 2927.5)/ 100000) )


resultArray = NULL
	
for(Dsim in 1:100) {
cat(Dsim, ".",sep="")
resultArray = abind(resultArray,
	PvalueMultiple(parameters, populationData, lambda, 
		c("Colon","Stomach"), 
  		verbose=F, deathFile=deathData, DiffCancer=DiffCancer), 
	along=6)
 }


pvalueSignif = resultArray[,,,,"pvalue",] < 0.05
# not sure if this will work
powers = apply(pvalueSignif, 1:4, mean, na.rm=T)


}


######################################
# p value function starts from here
#####################################
PvalueMultiple = function(parameters, populationData, lambda,
	CommonCancer, verbose=F, deathFile, DiffCancer) { 
# if oneEnvPerCommunity=T, each community has only one environment effect.  
#  otherwise there's one effect per person.									

	require(foreign)
	library(survival)
	

# simulate a cohort
datamat = simCohort(parameters, lambda, populationData, deathFile=deathFile, DiffCancer)


size=parameters$size
coefs = unlist(parameters[c("geno", "env", "geno:env")])

# create array for p values

Smethods = c("gaussian","gamma")

est=se = pvalue = array(1, c(length(Smethods), length(parameters$Followup), 
                    length(coefs), length(CommonCancer)),
		dimnames = list(Smethods, as.character(parameters$Followup), 
                      names(coefs), CommonCancer) )   

		
# if there's no one in the interaction groups, set p values to zero
if (sum(datamat$geno&datamat$env) == 0) {
	return(NULL)
} 


for(Dcancer in CommonCancer) {
	if(verbose) cat(Dcancer)
	thisCancer = rep(F, parameters$size)
# find the cancer we're interested in
  if (Dcancer=="All"){
  thisCancer = datamat$Cancer
  }else{                                  
	thisCancer[grep(Dcancer, datamat$CancerType)] = T
	}
	
	for(Dfollowup in parameters$Followup) {
	if(verbose) cat(Dfollowup)
    
	 	# set thisEventTime equal to age at end of followoup period
		ageEndFollowup = datamat$Age + Dfollowup - datamat$Enrollment
		hadEvent = datamat$Event < ageEndFollowup
		datamat$thisEventType = hadEvent &  thisCancer
		
		if(!any(datamat$thisEventType)) {
# if there is nobody with the interaction,  the pvalues to 1
  est[Smethod, as.character(Dfollowup),,Dcancer]= 	
    se[Smethod, as.character(Dfollowup),,Dcancer] = NA
		} else {	 # if there are  people with the interaction
 
		# censored failure time
		datamat$thisEventTime = datamat$Event
		datamat$thisEventTime[!hadEvent] = ageEndFollowup[!hadEvent]
# and add left truncation
    response = Surv(datamat$Age, datamat$thisEventTime, datamat$thisEventType)
	
		if(Dcancer == "Breast") {
warning("not impolimented")
		} else if (Dcancer=="Prostate") {
warning("not impolimented")
		}  

coxfitAfterYears =list()
		coxfitAfterYears$gamma=coxph(response~geno*env+
			strata(Gender,na.group=TRUE) + 
        frailty(CensusDivision, sparse=F), 
  		data=datamat, control=coxph.control(iter.max=50))

coxfitAfterYears$gaussian=coxph(response~geno*env+
			strata(Gender,na.group=TRUE) + 
      frailty(CensusDivision, distribution="gaussian", sparse=F), 
  		data=datamat, control=coxph.control(iter.max=50))
	
# p value is probability of z being above the observed z scores
# 1- probability of z being below the observed
for(Dmethod in Smethods)   { #loop through models
  est[Dmethod, as.character(Dfollowup),,Dcancer] <- 
	coxfitAfterYears[[Dmethod]]$coefficients[dimnames(est)[[3]] ]
allSe=	sqrt(diag(coxfitAfterYears[[Dmethod]]$var))
names(allSe) = names(coxfitAfterYears[[Dmethod]]$coefficients)
   se[Dmethod, as.character(Dfollowup),,Dcancer] <- allSe[dimnames(est)[[3]] ]

#pvalVec= 1 - pchisq((coef/ se)^2, 1)

#  pvalue[Dmethod, as.character(Dfollowup),,Dcancer]= 
#	pvalVec[dimnames(pvalue)[[3]] ]
	 	

} # end loop through models

	if(verbose) cat("\n")

} # end if have people with interaction

} # end followup years

} # end cancers 

# replace NA's with 1?
pvalue = 1 - pchisq((est/ se)^2, 1)


return(abind(pvalue = pvalue, coef=est, se=se, along = length(dim(pvalue))+1))
}





